{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/W2D1-postcourse-bugfix/tutorials/W2D2_LinearSystems/student/W2D2_Tutorial4.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "# Neuromatch Academy 2020, Week 2, Day 2, Tutorial 4\n",
    "\n",
    "# Autoregressive models\n",
    "\n",
    "**Content Creators**: Bing Wen Brunton, Biraj Pandey\n",
    "\n",
    "**Content Reviewers**: Norma Kuhn, John Butler, Matthew Krause, Ella Batty, Richard Gao, Michael Waskom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Tutorial Objectives\n",
    "\n",
    "The goal of this tutorial is to use the modeling tools and intuitions developed in the previous few tutorials and use them to _fit data_. The concept is to flip the previous tutorial -- instead of generating synthetic data points from a known underlying process, what if we are given data points measured in time and have to learn the underlying process?\n",
    "\n",
    "This tutorial is in two sections. \n",
    "\n",
    "**Section 1** walks through using regression of data to solve for the coefficient of an OU process from Tutorial 3. Next, **Section 2** generalizes this auto-regression framework to high-order autoregressive models, and we will try to fit data from monkeys at typewriters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "#@title Figure settings\n",
    "import ipywidgets as widgets       # interactive display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "# @title Helper Functions\n",
    "# drift-diffusion model, from Tutorial 3\n",
    "def ddm(T, x0, xinfty, lam, sig):\n",
    "  '''\n",
    "  Samples a trajectory of a drift-diffusion model.\n",
    "\n",
    "  args:\n",
    "  T (integer): length of time of the trajectory\n",
    "  x0 (float): position at time 0\n",
    "  xinfty (float): equilibrium position\n",
    "  lam (float): process param\n",
    "  sig: standard deviation of the normal distribution\n",
    "\n",
    "  returns:\n",
    "  t (numpy array of floats): time steps from 0 to T sampled every 1 unit\n",
    "  x (numpy array of floats): position at every time step\n",
    "  '''\n",
    "  t = np.arange(0, T, 1.)\n",
    "  x = np.zeros_like(t)\n",
    "  x[0] = x0\n",
    "\n",
    "  for k in range(len(t)-1):\n",
    "      x[k+1] = xinfty + lam * (x[k] - xinfty) + sig * np.random.standard_normal(size=1)\n",
    "\n",
    "  return t, x\n",
    "\n",
    "def build_time_delay_matrices(x, r):\n",
    "    \"\"\"\n",
    "    Builds x1 and x2 for regression\n",
    "\n",
    "    Args:\n",
    "    x (numpy array of floats): data to be auto regressed\n",
    "    r (scalar): order of Autoregression model\n",
    "\n",
    "    Returns:\n",
    "    (numpy array of floats) : to predict \"x2\"\n",
    "    (numpy array of floats) : predictors of size [r,n-r], \"x1\"\n",
    "\n",
    "    \"\"\"\n",
    "    # construct the time-delayed data matrices for order-r AR model\n",
    "    x1 = np.ones(len(x)-r)\n",
    "    x1 = np.vstack((x1, x[0:-r]))\n",
    "    xprime = x\n",
    "    for i in range(r-1):\n",
    "        xprime = np.roll(xprime, -1)\n",
    "        x1 = np.vstack((x1, xprime[0:-r]))\n",
    "\n",
    "    x2 = x[r:]\n",
    "\n",
    "    return x1, x2\n",
    "\n",
    "def AR_model(x, r):\n",
    "    \"\"\"\n",
    "    Solves Autoregression problem of order (r) for x\n",
    "\n",
    "    Args:\n",
    "    x (numpy array of floats): data to be auto regressed\n",
    "    r (scalar): order of Autoregression model\n",
    "\n",
    "    Returns:\n",
    "    (numpy array of floats) : to predict \"x2\"\n",
    "    (numpy array of floats) : predictors of size [r,n-r], \"x1\"\n",
    "    (numpy array of floats): coefficients of length [r] for prediction after\n",
    "    solving the regression problem \"p\"\n",
    "    \"\"\"\n",
    "    x1, x2 = build_time_delay_matrices(x, r)\n",
    "\n",
    "    # solve for an estimate of lambda as a linear regression problem\n",
    "    p, res, rnk, s = np.linalg.lstsq(x1.T, x2, rcond=None)\n",
    "\n",
    "    return x1, x2, p\n",
    "\n",
    "def AR_prediction(x_test, p):\n",
    "    \"\"\"\n",
    "    Returns the prediction for test data \"x_test\" with the regression\n",
    "    coefficients p\n",
    "\n",
    "    Args:\n",
    "    x_test (numpy array of floats): test data to be predicted\n",
    "    p (numpy array of floats): regression coefficients of size [r] after\n",
    "    solving the autoregression (order r) problem on train data\n",
    "\n",
    "    Returns:\n",
    "    (numpy array of floats): Predictions for test data. +1 if positive and -1\n",
    "    if negative.\n",
    "    \"\"\"\n",
    "    x1, x2 = build_time_delay_matrices(x_test, len(p)-1)\n",
    "\n",
    "    # Evaluating the AR_model function fit returns a number.\n",
    "    # We take the sign (- or +) of this number as the model's guess.\n",
    "    return np.sign(np.dot(x1.T, p))\n",
    "\n",
    "def error_rate(x_test, p):\n",
    "    \"\"\"\n",
    "    Returns the error of the Autoregression model. Error is the number of\n",
    "    mismatched predictions divided by total number of test points.\n",
    "\n",
    "    Args:\n",
    "    x_test (numpy array of floats): data to be predicted\n",
    "    p (numpy array of floats): regression coefficients of size [r] after\n",
    "    solving the autoregression (order r) problem on train data\n",
    "\n",
    "    Returns:\n",
    "    (float): Error (percentage).\n",
    "    \"\"\"\n",
    "    x1, x2 = build_time_delay_matrices(x_test, len(p)-1)\n",
    "\n",
    "    return np.count_nonzero(x2 - AR_prediction(x_test, p)) / len(x2)\n",
    "\n",
    "def plot_residual_histogram(res):\n",
    "  \"\"\"Helper function for Exercise 4A\"\"\"\n",
    "  fig = plt.figure()\n",
    "\n",
    "  plt.hist(res)\n",
    "  plt.xlabel('error in linear model')\n",
    "  plt.title('stdev of errors = {std:.4f}'.format(std=res.std()))\n",
    "\n",
    "  plt.show()\n",
    "\n",
    "def plot_training_fit(x1, x2, p):\n",
    "  \"\"\"Helper function for Exercise 4B\"\"\"\n",
    "  fig = plt.figure()\n",
    "\n",
    "  plt.scatter(x2 + np.random.standard_normal(len(x2))*0.02,\n",
    "              np.dot(x1.T, p), alpha=0.2)\n",
    "  plt.title('Training fit, order {r:d} AR model'.format(r=r))\n",
    "\n",
    "  plt.xlabel('x')\n",
    "  plt.ylabel('estimated x')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "# Section 1: Fitting data to the OU process\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 519
    },
    "colab_type": "code",
    "outputId": "9918203f-5182-4350-8b44-319d14d1ea78"
   },
   "outputs": [],
   "source": [
    "#@title Video 1: Autoregressive models\n",
    "# Insert the ID of the corresponding youtube video\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"VdiVSTPbJ7I\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "To see how this works, let's continue the previous example with the drift-diffusion (OU) process. Our process had the following form:\n",
    "\n",
    "$x_{k+1} = x_{\\infty} + \\lambda(x_k - x_{\\infty}) + \\sigma \\eta$\n",
    "\n",
    "where $\\eta$ is sampled from a standard normal distribution. \n",
    "\n",
    "For simplity, we set $x_\\infty = 0$. Let's plot a trajectory for this process again below. Take note of the parameters of the process because they will be important later. **Run the code cell below.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "outputId": "cc15d0f7-6d54-43c4-e99e-7a4470d6f372"
   },
   "outputs": [],
   "source": [
    "#@title Simulating the drift diffusion model\n",
    "np.random.seed(2020) # set random seed\n",
    "\n",
    "# parameters\n",
    "T = 200\n",
    "x0 = 10\n",
    "xinfty = 0\n",
    "lam = 0.9\n",
    "sig = 0.2\n",
    "\n",
    "# drift-diffusion model from tutorial 3\n",
    "t, x = ddm(T, x0, xinfty, lam, sig)\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.title('$x_0=%d, x_{\\infty}=%d, \\lambda=%0.1f, \\sigma=%0.1f$' % (x0, xinfty, lam, sig))\n",
    "plt.plot(t, x, 'k.')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('position x')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "What if we were given these positions $x$ as they evolve in time as data, how would we get back out the dynamics of the system $\\lambda$? \n",
    "\n",
    "Since a little bird told us that this system takes on the form\n",
    "\n",
    "$x_{k+1} = \\lambda x_k + \\eta$,\n",
    "\n",
    "where $\\eta$ is noise from a normal distribution, our approach is to solve for $\\lambda$ as a **regression problem**. \n",
    "\n",
    "As a check, let's plot every pair of points adjacent in time ($x_{k+1}$ vs. $x_k$) against eachother to see if there is a linear relationship between them. **Run the code cell below.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 431
    },
    "colab_type": "code",
    "outputId": "f01962a8-e4c5-471c-b03d-0fcb194ec78a"
   },
   "outputs": [],
   "source": [
    "# @title X(k) vs. X(k+1)\n",
    "# make a scatter plot of every data point in x\n",
    "# at time k versus time k+1\n",
    "fig = plt.figure()\n",
    "plt.scatter(x[0:-2], x[1:-1], color='k')\n",
    "plt.plot([0, 10], [0, 10], 'k--', label='$x_{k+1} = x_k$ line')\n",
    "plt.xlabel('$x_k$')\n",
    "plt.ylabel('$x_{k+1}$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Hooray, it's a line! This is evidence that that the _dynamics that generated the data_ is **linear**. We can now reformulate this task as a regression problem.\n",
    "\n",
    "Let $\\mathbf{x_1} = x_{0:T-1}$ and $\\mathbf{x_2} = x_{1:T}$ be vectors of the data indexed so that they are shifted in time by one. Then, our regression problem is\n",
    "\n",
    "$$\\mathbf{x}_2 = \\lambda \\mathbf{x}_1$$\n",
    "\n",
    "This model is **autoregressive**, where _auto_ means self. In other words, it's a regression of the time series on itself from the past. The equation as written above is only a function of itself from _one step_ in the past, so we can call it a _first order_ autoregressive model.\n",
    "\n",
    "Now, let's set up the regression problem below and solve for $\\lambda.$ We will plot our data with the regression line to see if they agree. **Run the code cell below.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "outputId": "a972266a-b33e-4340-b904-00b3ff03c896"
   },
   "outputs": [],
   "source": [
    "#@title Solving for lambda through autoregression\n",
    "# build the two data vectors from x\n",
    "x1 = x[0:-2]\n",
    "x1 = x1[:, np.newaxis]**[0, 1]\n",
    "\n",
    "x2 = x[1:-1]\n",
    "\n",
    "# solve for an estimate of lambda as a linear regression problem\n",
    "p, res, rnk, s = np.linalg.lstsq(x1, x2, rcond=None)\n",
    "\n",
    "# here we've artificially added a vector of 1's to the x1 array,\n",
    "# so that our linear regression problem has an intercept term to fit.\n",
    "# we expect this coefficient to be close to 0.\n",
    "# the second coefficient in the regression is the linear term:\n",
    "# that's the one we're after!\n",
    "lam_hat = p[1]\n",
    "\n",
    "# plot the data points\n",
    "fig = plt.figure()\n",
    "plt.scatter(x[0:-2], x[1:-1], color='k')\n",
    "plt.xlabel('$x_k$')\n",
    "plt.ylabel('$x_{k+1}$')\n",
    "\n",
    "# plot the 45 degree line\n",
    "plt.plot([0, 10], [0, 10], 'k--', label='$x_{k+1} = x_k$ line')\n",
    "\n",
    "\n",
    "# plot the regression line on top\n",
    "xx = np.linspace(-sig*10, max(x), 100)\n",
    "yy = p[0] + lam_hat * xx\n",
    "plt.plot(xx, yy, 'r', linewidth=2, label='regression line')\n",
    "\n",
    "mytitle = 'True $\\lambda$ = {lam:.4f}, Estimate $\\lambda$ = {lam_hat:.4f}'\n",
    "plt.title(mytitle.format(lam=lam, lam_hat=lam_hat))\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Pretty cool! So now we have a way to predict $x_{k+1}$ if given any data point $x_k$. Let's take a look at how accurate this one-step prediction might be by plotting the residuals.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 1 (4A): Residuals of the autoregressive model\n",
    "\n",
    "Plot a histogram of residuals of our autoregressive model, by taking the difference between the _data_ $\\mathbf{x_2}$ and the _model_ prediction. Do you notice anything about the standard deviation of these residuals and the equations that generated this synthetic dataset?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "code",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "##############################################################################\n",
    "## Insert your code here take to compute the residual (error)\n",
    "##############################################################################\n",
    "# compute the predicted values using the autoregressive model (lam_hat), and\n",
    "# the residual is the difference between x2 and the prediction\n",
    "# res = ...\n",
    "\n",
    "# Uncomment once you fill out above\n",
    "#plot_residual_histogram(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 433
    },
    "colab_type": "text",
    "outputId": "1e8c92c4-d29d-49a1-c023-d7400ee23437"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D2_LinearSystems/solutions/W2D2_Tutorial4_Solution_cd2c0633.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=416 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W2D2_LinearSystems/static/W2D2_Tutorial4_Solution_cd2c0633_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 2: Higher order autoregressive models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 519
    },
    "colab_type": "code",
    "outputId": "4066b8a2-d1f8-41ef-8c36-1cd30389340c"
   },
   "outputs": [],
   "source": [
    "#@title Video 2: Monkey at a typewriter\n",
    "# Insert the ID of the corresponding youtube video\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"f2z0eopWB8Y\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Now that we have established the autoregressive framework, generalizing for dependence on data points from the past is straightfoward. **Higher order** autoregression models a future time point based on _more than one points in the past_.\n",
    "\n",
    "In one dimension, we can write such an order-$r$ model as\n",
    "\n",
    "$x_{k+1} = \\alpha_0 + \\alpha_1 x_k + \\alpha_2 x_{k-1} + \\alpha_3 x_{k-2} + \\dots + \\alpha_{r+1} x_{k-r}$,\n",
    "\n",
    "where the $\\alpha$'s are the $r+1$ coefficients to be fit to the data available."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "These models are useful to account for some **history dependence** in the trajectory of timeseries. This next part of the tutorial will explore one such timeseries, and you can do an experiment on yourself!\n",
    "\n",
    "In particular, we will explore a binary random sequence of 0's and 1's that would occur if you flipped a coin and jotted down the flips. \n",
    "\n",
    "The difference is that, instead of actually flipping a coin (or using code to generate such a sequence), you -- yes you, human -- are going to generate such a random Bernoulli sequence as best as you can by typing in 0's and 1's. We will then build higher-order AR models to see if we can identify predictable patterns in the time-history of digits you generate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "**But first**, let's try this on a sequence with a simple pattern, just to make sure the framework is functional. Below, we generate an entirely predictable sequence and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "outputId": "b260e7b7-f0cf-440a-b572-7e8d083ba1d4"
   },
   "outputs": [],
   "source": [
    "# this sequence is entirely predictable, so an AR model should work\n",
    "monkey_at_typewriter = '1010101010101010101010101010101010101010101010101'\n",
    "\n",
    "# Bonus: this sequence is also predictable, but does an order-1 AR model work?\n",
    "#monkey_at_typewriter = '100100100100100100100100100100100100100'\n",
    "\n",
    "# function to turn chars to numpy array,\n",
    "# coding it this way makes the math easier\n",
    "# '0' -> -1\n",
    "# '1' -> +1\n",
    "def char2array(s):\n",
    "    m = [int(c) for c in s]\n",
    "    x = np.array(m)\n",
    "    return x*2 - 1\n",
    "\n",
    "x = char2array(monkey_at_typewriter)\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.step(x, '.-')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('random variable')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Now, let's set up our regression problem (order 1 autoregression like above) by defining $\\mathbf{x_1}$ and $\\mathbf{x_2}$ and solve it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "# build the two data vectors from x\n",
    "x1 = x[0:-2]\n",
    "x1 = x1[:, np.newaxis]**[0, 1]\n",
    "\n",
    "x2 = x[1:-1]\n",
    "\n",
    "# solve for an estimate of lambda as a linear regression problem\n",
    "p, res, rnk, s = np.linalg.lstsq(x1, x2, rcond=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "outputId": "337653d4-5874-43f7-99a2-992ef7e79cc8"
   },
   "outputs": [],
   "source": [
    "# take a look at the resulting regression coefficients\n",
    "print('alpha_0 = {a0:.2f}, alpha_1 = {a1:.2f}'.format(a0=p[0], a1=p[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Think:\n",
    "Do the values we got for $\\alpha_0$ and $\\alpha_1$ make sense? Write down the corresponding autoregressive model and convince yourself that it gives the alternating 0's and 1's we asked it to fit as data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D2_LinearSystems/solutions/W2D2_Tutorial4_Solution_3568b621.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Truly random sequences of numbers have no structure and should not be predictable by an AR or any other models.\n",
    "\n",
    "However, humans are notoriously terrible at generating random sequences of numbers! (Other animals are no better...)\n",
    "\n",
    "To test out an application of higher-order AR models, let's use them to **model a sequence of 0's and 1's that a human tried to produce at random**. In particular, I convinced my 9-yr-old monkey to sit at a typewriter (my laptop) and enter some digits as randomly as he is able. The digits he typed in are in the code, and we can plot them as a timeseries of digits here.\n",
    "\n",
    "If the digits really have no structure, then we expect our model to do about as well as guessing, producing an error rate of 0.5. Let's see how well we can do!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "outputId": "82f78b1c-b086-4f39-bfec-956e6848964a"
   },
   "outputs": [],
   "source": [
    "# data generated by 9-yr-ld JAB:\n",
    "# we will be using this sequence to train the data\n",
    "monkey_at_typewriter = '10010101001101000111001010110001100101000101101001010010101010001101101001101000011110100011011010010011001101000011101001110000011111011101000011110000111101001010101000111100000011111000001010100110101001011010010100101101000110010001100011100011100011100010110010111000101'\n",
    "\n",
    "# we will be using this sequence to test the data\n",
    "test_monkey = '00100101100001101001100111100101011100101011101001010101000010110101001010100011110'\n",
    "\n",
    "x = char2array(monkey_at_typewriter)\n",
    "test = char2array(test_monkey)\n",
    "\n",
    "## testing: machine generated randint should be entirely unpredictable\n",
    "## uncomment the lines below to try random numbers instead\n",
    "# np.random.seed(2020) # set random seed\n",
    "# x = char2array(np.random.randint(2, size=500))\n",
    "# test = char2array(np.random.randint(2, size=500))\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.step(x, '.-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 2 (4B): Fitting AR models\n",
    "\n",
    "Fit a order-5 ($r=5$) AR model to the data vector $x$. To do this, we have included some helper functions, including ``AR_model``. \n",
    "\n",
    "We will then plot the observations against the trained model. Note that this means we are using a sequence of the previous 5 digits to predict the next one. \n",
    "\n",
    "Additionally, output from our regression model are continuous (real numbers) whereas our data are scalar (+1/-1). So, we will take the sign of our continuous outputs (+1 if positive and -1 if negative) as our predictions to make them comparable with data. Our error rate will simply be the number of mismatched predictions divided by the total number of predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 278
    },
    "colab_type": "code",
    "outputId": "7dd3e4df-7f91-4e8b-f9ad-2fc70f5daa04"
   },
   "outputs": [],
   "source": [
    "# Let's see what our function AR model entails\n",
    "help(AR_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "code",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "##############################################################################\n",
    "## TODO: Insert your code here for fitting the AR model\n",
    "##############################################################################\n",
    "# define the model order, and use AR_model() to generate the model and prediction\n",
    "# r = ...\n",
    "# x1, x2, p = AR_model(...)\n",
    "\n",
    "\n",
    "# Uncomment below once you've completed above\n",
    "# Plot the Training data fit\n",
    "# Note that this adds a small amount of jttter to horizontal axis for visualization purposes\n",
    "# plot_training_fit(x1, x2, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 433
    },
    "colab_type": "text",
    "outputId": "0d6af7b2-0a87-435f-c0b5-a55bc7ec8c72"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D2_LinearSystems/solutions/W2D2_Tutorial4_Solution_7171e02c.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=416 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W2D2_LinearSystems/static/W2D2_Tutorial4_Solution_7171e02c_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Let's check out how the model does on the test data that it's never seen before!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 447
    },
    "colab_type": "code",
    "outputId": "1c0b8fac-1b50-48f4-a97b-329c0fd59ff6"
   },
   "outputs": [],
   "source": [
    "x1_test, x2_test = build_time_delay_matrices(test, r)\n",
    "fig = plt.figure()\n",
    "plt.scatter(x2_test+np.random.standard_normal(len(x2_test))*0.02,\n",
    "            np.dot(x1_test.T, p), alpha=0.5)\n",
    "\n",
    "mytitle = 'Testing fit, order {r:d} AR model, err = {err:.3f}'\n",
    "plt.title(mytitle.format(r=r, err=error_rate(test, p)))\n",
    "\n",
    "plt.xlabel('test x')\n",
    "plt.ylabel('estimated x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Not bad! We're getting errors that are smaller than 0.5 (what we would have gotten by chance).\n",
    "\n",
    "Let's now try **AR models of different orders** systematically, and plot the test error of each.\n",
    "\n",
    "_Remember_: The model has never seen the test data before, and random guessing would produce an error of $0.5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "outputId": "cfa6c0e3-687e-494c-a58e-8b9f9131d1b1"
   },
   "outputs": [],
   "source": [
    "# range of r's to try\n",
    "r = np.arange(1, 21)\n",
    "err = np.ones_like(r) * 1.0\n",
    "\n",
    "for i, rr in enumerate(r):\n",
    "    # fitting the model on training data\n",
    "    x1, x2, p = AR_model(x, rr)\n",
    "    # computing and storing the test error\n",
    "    test_error = error_rate(test, p)\n",
    "    err[i] = test_error\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.plot(r, err, '.-')\n",
    "plt.plot([1, r[-1]], [0.5, 0.5], c='r', label='random chance')\n",
    "plt.xlabel('Order r of AR model')\n",
    "plt.ylabel('Test error')\n",
    "plt.xticks(np.arange(0,25,5))\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Notice that there's a sweet spot in the test error! The 6th order AR model does a really good job here, and for larger $r$'s, the model starts to overfit the training data and does not do well on the test data.\n",
    "\n",
    "In summary:\n",
    "\n",
    "\"**I can't believe I'm so predictable!**\" - JAB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "In this tutorial, we learned:\n",
    "\n",
    "* How learning the parameters of a linear dynamical system can be formulated as a regression problem from data.\n",
    "* Time-history dependence can be incorporated into the regression framework as a multiple regression problem.\n",
    "* That humans are no good at generating random (not predictable) sequences. Try it on yourself!"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W2D2_Tutorial4",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
